{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run 4.3.2-common.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Principal values are tensor([9736.4043, 7876.6597])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAHHCAYAAAC1G/yyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOzde1zUdb4/8NcwCiIyIMpNQfFSiml5tNa0KNhINPToIrpaW+jaZV1tQa3W2vLSjdI6QW5p9TuPdGv1pC7ZbmsWFZR7ZK1My7t1glAUxBsjqKAz398f43ecgbl8B2a+t3k9e/hIvvNl5jPgzPc1n8/78/kYBEEQQEREREQAgBClG0BERESkJgxHRERERA4YjoiIiIgcMBwREREROWA4IiIiInLAcERERETkgOGIiIiIyAHDEREREZEDhiMiIiIiBwxHRKRpBoMBS5cuVboZTmbOnImUlJR2f39KSgpmzpzpt/YQkW8Yjog0aunSpTAYDDh58qTL24cOHYr09HR5G6Uh+/fvx9KlS1FVVaV0U/zm2LFjWLp0KXbv3q10UwAAW7ZsUV1wJZKC4YiIgtL+/fuxbNky3YWjZcuWqSocLVu2TOlmEPmM4YiIZNfU1OTyuCAIuHDhgsytISJyxnBEFCTKy8thMBiwYcMGPPfcc0hKSkKXLl1wxx134Mcff2xz/o4dO3DXXXehe/fuiIiIwPXXX4/i4mKncz7//HOkpaUhIiIC0dHRmDRpEg4cOOB0jjj8t3//ftx9993o3r07br31VgC22poJEybg448/xo033ojw8HC88cYbAICzZ8+ioKAAycnJCAsLw8CBA/Hiiy/CarV6fJ4///wzfv/732PQoEEIDw9Hjx49MHXqVKceojVr1mDq1KkAgIyMDBgMBhgMBpSXl9vP+eijj+zPLTIyEtnZ2di3b1+bx9u8eTOGDh2KLl26YOjQoXj//fc9ts+RIAh49tlnkZSUhK5duyIjI8PlY5w+fRqPPPIIhg0bhm7dusFkMmH8+PH47rvv7OeUl5fjpptuAgDMmjXL/pzWrFkDANi2bRumTp2KPn36ICwsDMnJyZg/f36bMFpbW4tZs2YhKSkJYWFhSExMxKRJk9r0sHn7+cycOROvvfYaANjbYjAYJP9siJTUSekGEJG8XnjhBYSEhOCRRx5BQ0MDli9fjnvuuQc7duywn1NaWooJEyYgMTER+fn5SEhIwIEDB/Dhhx8iPz8fAPDpp59i/Pjx6N+/P5YuXYoLFy5g5cqVuOWWW/Dtt9+2KUieOnUqrrnmGjz//PMQBMF+/NChQ5gxYwYeeughPPDAAxg0aBDOnz+P22+/HTU1NXjooYfQp08fbN++HY8//jiOHz+OoqIit8/v66+/xvbt2zF9+nQkJSWhqqoKq1atQnp6Ovbv34+uXbvitttuwx/+8Ae8+uqreOKJJ5CamgoA9v+/8847yMvLQ1ZWFl588UWcP38eq1atwq233opdu3bZn9snn3yCKVOmYMiQISgsLMSpU6fswUKKxYsX49lnn8Vdd92Fu+66C99++y3Gjh2LlpYWp/N++uknbN68GVOnTkW/fv1QV1eHN954A7fffjv279+PXr16ITU1FU8//TQWL16MBx98EGlpaQCAMWPGAAA2btyI8+fPY86cOejRowe++uorrFy5EkePHsXGjRvtjzVlyhTs27cPDz/8MFJSUnDixAmUlpaiurra/ryl/HweeughHDt2DKWlpXjnnXck/TyIVEMgIk1asmSJAECor693eft1110n3H777favy8rKBABCamqq0NzcbD9eXFwsABD27NkjCIIgXL58WejXr5/Qt29f4cyZM073abVa7X8fPny4EBcXJ5w6dcp+7LvvvhNCQkKE++67r007Z8yY0aaNffv2FQAIW7dudTr+zDPPCBEREcLhw4edji9atEgwGo1CdXW1/RgAYcmSJfavz58/3+ZxKioqBADCX/7yF/uxjRs3CgCEsrIyp3PPnTsnREdHCw888IDT8draWiEqKsrp+PDhw4XExETh7Nmz9mOffPKJAEDo27dvm3Y4OnHihBAaGipkZ2c7/VyfeOIJAYCQl5dnP3bx4kXBYrE4fX9lZaUQFhYmPP300/ZjX3/9tQBAePvtt9s8nqufS2FhoWAwGISff/5ZEARBOHPmjABAWLFihdt2+/LzmTt3rsDLDGkRh9WIgsysWbMQGhpq/1rsYfjpp58AALt27UJlZSUKCgoQHR3t9L3isMjx48exe/duzJw5EzExMfbbr7/+etx5553YsmVLm8f93e9+57I9/fr1Q1ZWltOxjRs3Ii0tDd27d8fJkyftfzIzM2GxWPDll1+6fX7h4eH2v1+6dAmnTp3CwIEDER0djW+//dbt94lKS0tx9uxZzJgxw+mxjUYjRo0ahbKyMqefQV5eHqKiouzff+edd2LIkCFeH+fTTz9FS0sLHn74YafhpoKCgjbnhoWFISTE9nZtsVhw6tQpdOvWDYMGDZL0nADnn0tTUxNOnjyJMWPGQBAE7Nq1y35OaGgoysvLcebMGZf3I/XnQ6RlHFYj0jFXNR59+vRx+rp79+4AYL8Y/t///R8A21IA7vz8888AgEGDBrW5LTU1FR9//DGampoQERFhP96vXz+X9+Xq+A8//IDvv/8esbGxLr/nxIkTbtt24cIFFBYW4u2330ZNTY3TEF5DQ4Pb73N8bAD45S9/6fJ2k8kE4OrP4JprrmlzjpTQ4u77Y2Nj7b8TkdVqRXFxMV5//XVUVlbCYrHYb+vRo4fHxxFVV1dj8eLF+Pvf/94m+Ig/l7CwMLz44otYuHAh4uPjcfPNN2PChAm47777kJCQAED6z4dIyxiOiDSqS5cuAOB2dtf58+ft5zgyGo0uz3cMEYHg2HPh7bjVasWdd96Jxx57zOX3XHvttW4f5+GHH8bbb7+NgoICjB49GlFRUTAYDJg+fbrXYm7xsQFbXY0YCBx16iT/2+bzzz+Pp556Cr/97W/xzDPPICYmBiEhISgoKJD0nCwWC+68806cPn0af/zjHzF48GBERESgpqYGM2fOdLqPgoICTJw4EZs3b8bHH3+Mp556CoWFhfj888/xH//xH6r8+RD5G/8VE2lU3759AdgKmpOTk51uO3/+PI4cOYKxY8f6fL8DBgwAAOzduxeZmZleH7u1gwcPomfPnk69Ru1pQ2Njo9vH92TTpk3Iy8vDyy+/bD928eJFnD171uk8dzOnxOcfFxfn8fHFn4HYk+LI1c/F0/f379/ffry+vr5Nz86mTZuQkZGB//7v/3Y6fvbsWfTs2dP+tbvntGfPHhw+fBhr167FfffdZz9eWlrq8vwBAwZg4cKFWLhwIX744QcMHz4cL7/8Mt59913JPx9P7SFSO9YcEWnUHXfcgdDQUKxatapN78Gbb76Jy5cvY/z48T7f74gRI9CvXz8UFRW1CRRi71JiYiKGDx+OtWvXOp2zd+9efPLJJ7jrrrva8YyumjZtGioqKvDxxx+3ue3s2bO4fPmy2+81Go1tesFWrlzpNBQFwB7eWj/HrKwsmEwmPP/887h06VKb+6+vrwfg/DNwHK4rLS3F/v37vTxDIDMzE507d8bKlSud2utqJp6r57Rx40bU1NRIek5ib6HjfQiC0GZphvPnz+PixYtOxwYMGIDIyEg0NzcDkP7z8dQeIrVjzxGRRsXFxWHx4sV48skncdttt+E///M/0bVrV2zfvh3r16/H2LFjMXHiRJ/vNyQkBKtWrcLEiRMxfPhwzJo1C4mJiTh48CD27dtnDywrVqzA+PHjMXr0aMyePds+lT8qKqrDW0Y8+uij+Pvf/44JEyZg5syZGDlyJJqamrBnzx5s2rQJVVVVTj0mjiZMmIB33nkHUVFRGDJkCCoqKvDpp5+2qc0ZPnw4jEYjXnzxRTQ0NCAsLAy//OUvERcXh1WrVuHee+/FiBEjMH36dMTGxqK6uhr//Oc/ccstt+DPf/4zAKCwsBDZ2dm49dZb8dvf/hanT5/GypUrcd1116GxsdHjc4yNjcUjjzyCwsJCTJgwAXfddRd27dqFjz76qM1zmzBhAp5++mnMmjULY8aMwZ49e/DXv/7VqccJsAWZ6OhorF69GpGRkYiIiMCoUaMwePBgDBgwAI888ghqampgMpnwt7/9rU0P1eHDh3HHHXdg2rRpGDJkCDp16oT3338fdXV1mD59OgBbTZHUn8/IkSMBAH/4wx+QlZUFo9Fovx8iVVNsnhwR+cW7774r3HzzzUJERIQQFhYmDB48WFi2bJlw8eJFp/PEqfwbN250Ol5ZWely+ve//vUv4c477xQiIyOFiIgI4frrrxdWrlzpdM6nn34q3HLLLUJ4eLhgMpmEiRMnCvv373c6x9OSA3379hWys7NdPq9z584Jjz/+uDBw4EAhNDRU6NmzpzBmzBjhpZdeElpaWuznodVU/jNnzgizZs0SevbsKXTr1k3IysoSDh48KPTt29dperwgCMJbb70l9O/fXzAajW2m9ZeVlQlZWVlCVFSU0KVLF2HAgAHCzJkzhW+++cbpPv72t78JqampQlhYmDBkyBChpKREyMvL8zqVXxAEwWKxCMuWLRMSExOF8PBwIT09Xdi7d2+btl68eFFYuHCh/bxbbrlFqKioEG6//Xan5RoEQRA++OADYciQIUKnTp2cfq/79+8XMjMzhW7dugk9e/YUHnjgAeG7775zOufkyZPC3LlzhcGDBwsRERFCVFSUMGrUKGHDhg1t2i7l53P58mXh4YcfFmJjYwWDwcBp/aQZBkEIcBUmERERkYaw5oiIiIjIAcMRERERkQOGIyIiIiIHDEdEREREDhiOiIiIiBwwHBERERE54CKQPrJarTh27BgiIyO5ND4REZFGCIKAc+fOoVevXggJ8dw3xHDko2PHjrXZx4qIiIi04ciRI0hKSvJ4DsORjyIjIwHYfrgmk0nh1hAREZEUZrMZycnJ9uu4JwxHPhKH0kwmE8MRERGRxkgpiWFBNhEREZEDhiMiIiIiBwxHRERERA4YjoiIiIgcMBwREREROWA4IiIiInLAcERERETkgOGIiIiIyAHDEREREZEDrpBNREQuWawWbKvehuPnjiMxMhFpfdJgDDEq3SyigGM4IiKiNkoOlCB/az6Omo/ajyWZklA8rhg5qTkKtowo8DisRkQURCxWC8qryrF+z3qUV5XDYrW0OafkQAlyN+Q6BSMAqDHXIHdDLkoOlMjVXCJFaCYcFRYW4qabbkJkZCTi4uIwefJkHDp0yOmcixcvYu7cuejRowe6deuGKVOmoK6uzumc6upqZGdno2vXroiLi8Ojjz6Ky5cvy/lUiIgUUXKgBCnFKchYm4G7S+5GxtoMpBSnOIUdi9WC/K35ECC0+X7xWMHWApehikgvNBOOvvjiC8ydOxf//ve/UVpaikuXLmHs2LFoamqynzN//nz84x//wMaNG/HFF1/g2LFjyMm52v1rsViQnZ2NlpYWbN++HWvXrsWaNWuwePFiJZ4SEZFspPYGbave1uYcRwIEHDEfwbbqbQFtL5GSDIIgtP14oAH19fWIi4vDF198gdtuuw0NDQ2IjY3FunXrkJubCwA4ePAgUlNTUVFRgZtvvhkfffQRJkyYgGPHjiE+Ph4AsHr1avzxj39EfX09QkNDvT6u2WxGVFQUGhoaYDKZAvociYj8wWK1IKU4xW3oMcCAJFMSKvMrsWHfBtxdcrfX+1yXsw4zhs3wd1OJAsaX67dmeo5aa2hoAADExMQAAHbu3IlLly4hMzPTfs7gwYPRp08fVFRUAAAqKiowbNgwezACgKysLJjNZuzbt8/l4zQ3N8NsNjv9ISLSEl96gxIjEyXdp9TziLRIk+HIarWioKAAt9xyC4YOHQoAqK2tRWhoKKKjo53OjY+PR21trf0cx2Ak3i7e5kphYSGioqLsf5KTk/39dIiIAur4ueOSz0vrk4YkUxIMMLg8xwADkk3JSOuT5s8mEqmKJsPR3LlzsXfvXvzP//xPwB/r8ccfR0NDg/3PkSNHAv6YRET+5EtvkDHEiOJxxQDQJiCJXxeNK+J6R6RrmgtH8+bNw4cffoiysjIkJSXZjyckJKClpQVnz551Or+urg4JCQn2c1rPXhO/Fs9pLSwsDCaTyekPEZGW+NoblJOag03TNqG3qbfTeUmmJGyatonrHJHuaSYcCYKAefPm4f3338fnn3+Ofv36Od0+cuRIdO7cGZ999pn92KFDh1BdXY3Ro0cDAEaPHo09e/bgxIkT9nNKS0thMpkwZMgQeZ4IEZHM2tMblJOag6r8KpTllWFdzjqU5ZWhMr+SwYiCgmZmq/3+97/HunXr8MEHH2DQoEH241FRUQgPDwcAzJkzB1u2bMGaNWtgMpnw8MMPAwC2b98OwDaVf/jw4ejVqxeWL1+O2tpa3Hvvvbj//vvx/PPPS2oHZ6sRkVa5WvU62ZSMonFFDD2ke75cvzUTjgwG193Bb7/9NmbOnAnAtgjkwoULsX79ejQ3NyMrKwuvv/6605DZzz//jDlz5qC8vBwRERHIy8vDCy+8gE6dpO2kwnBERFrG/dIoWOkyHKkFwxEREZH2BMU6R0RERESBwHBERERE5IDhiIiIiMgBwxERERGRA4YjIiIiIgcMR0REREQOGI6IiIiIHDAcERERETlgOCIiIiJywHBERERE5IDhiIiIiMgBwxERERGRA4YjIiIiIgcMR0REREQOGI6IiIiIHHRSugFERB1lsVqwrXobjp87jsTIRKT1SYMxxKh0s4hIoxiOiEjTSg6UIH9rPo6aj9qPJZmSUDyuGDmpOQq2jIi0isNqRKRZJQdKkLsh1ykYAUCNuQa5G3JRcqBEoZYRkZYxHBGRJlmsFuRvzYcAoc1t4rGCrQWwWC1yN42INI7hiIg0aVv1tjY9Ro4ECDhiPoJt1dtkbBUR6QHDERFp0vFzx/16HhGRiOGIiDQpMTLRr+cREYkYjohIk9L6pCHJlAQDDC5vN8CAZFMy0vqkydwyItI6hiMi0qwHRjzgsiBbDExF44q43hER+YzrHBGR5rha28hRkikJReOKuM4REbULwxERaYq4tpGrHiMAWJa+DH9K+xN7jIio3TisRkSa4WltI8A2nPb/vv1/MreKiPSG4YiININrGxGRHBiOiEgzuLYREcmB4YiINKPdaxstWgTk5ACHDwegVUSkNwxHRKQZ7V7baMsW4P33gaqqwDeSiDSP4YiINMMYYkTxuGIAaBOQPK5tFBdn+399fcDbSETax3BERJqSk5qDTdM2obept9PxJFMSNk3b5Hpto9hY2//r6mRoIRFpHdc5IiLNyUnNwaRBk7CtehuOnzuOxMhEpPVJc7+2EXuOiMgHDEdEpEnGECPSU9KlnSyGI/YcEZEEHFYjIv0Th9XYc0REErDniIhcslgt0oet1I7DakTkA4YjImrD1cauSaYkFI8r1uZmrhxWIyIfcFiNiJyIG7u23qajxlyD3A25KDlQolDLOoA9R0TkA4YjIrLztLGreKxgawEsVovcTesYseaoqQk4f17ZthCR6jEcEZGdbjd2NZmA0FDb39l75JbFakF5VTnW71mP8qpy7YVgIj9hzRER2el2Y1eDwTa0dvQocOIE0Lev0i1SHd3VmRF1AHuOiMiu3Ru7agHrjtzSZZ0ZUQcwHBGRXbs3dtUCse7oxAll26Eyuq0zI+oAhiMismv3xq5awJ4jl3RbZ0bUAQxHROSkXRu7aoEYjthz5ES3dWZEHcCCbCJqw+eNXbWAw2ou6brOjKidGI6IyCWfNnbVAg6ruSTWmdWYa1zWHRlgQJIpSZt1ZkTtxGE1IgoO3ELEJV3XmRG1E8MREQUH9hy5pds6M6J24rAaEQUHx5ojQbAtDEl2uqwzI2onhiMiCg5iOGpuBhobgchIZdujQrqrMyNqJ4YjIgoOERG2P01NtrojhiNNs1gt7OWigGE4IqLgERcHVFba6o4GDlS6NdRO3AeOAo0F2UQ6xl3WWxGH1liUrVncB47kwJ4jIh2yWC14bttzKN5RjNMXTtuPB/2na07n1zRv+8AZYEDB1gJMGjSJQ2zUIew5ItKZkgMliH8pHkvKlzgFI4CfrjmdX9u4DxzJheGISEdKDpRgyoYpOHXhlMvbg36Xde6vpmncB47kwnBEpBPikIM3Qf3pmjVHmsZ94EgurDki0glvQw6tBeWnaxXXHHFqunfcB47kwnBEpBO+hp2g/HSt0pojtU5NV1tgE/eBy92QCwMMTgHJ3/vAqe25k7wYjoh0wpewk2xK1s2na58uYleG1YQTJ/BFVbkqLnzi1PTWPSFi8bxSe5upNbCJ+8C5alvRuCK/tE2tz53kYxAEoW3fJLllNpsRFRWFhoYGmEwmpZtDZGexWpBSnOJ2yEFkgEE3m4n6fBGrqQGSknApBAh7EhBCJHxPAIm/M3fDoeIwUWV+pazhzV1gE3tn1PDvJ1A9O1p47tQ+vly/NVWQ/eWXX2LixIno1asXDAYDNm/e7HS7IAhYvHgxEhMTER4ejszMTPzwww9O55w+fRr33HMPTCYToqOjMXv2bDQ2Nsr5NIgCQhxyAK6+kbfWI7yHbt7c27MY4Psn/wUA6GwFoi9K+55AUuPUdG9rCQHqmO0o7gM3Y9gMpKek+20oTQvPnQJPU+GoqakJN9xwA1577TWXty9fvhyvvvoqVq9ejR07diAiIgJZWVm4ePHqu+A999yDffv2obS0FB9++CG+/PJLPPjgg3I9BaKAEoccept6Ox2PCY/BsvRlqHukThfByJeLmLhK+F+//yse/GQezobZzotrcv89clHj1HQ1Bja5BPNzJ2eaqjkaP348xo8f7/I2QRBQVFSEJ598EpMmTQIA/OUvf0F8fDw2b96M6dOn48CBA9i6dSu+/vpr3HjjjQCAlStX4q677sJLL72EXr16yfZciAIlJzUHkwZNkq2YVInCVakXsee2PYe3vn3L6dwTEUB0MxB7Hjjk4nu2VW+TbWd6NU5NV2Ngk0swP3dypqlw5EllZSVqa2uRmZlpPxYVFYVRo0ahoqIC06dPR0VFBaKjo+3BCAAyMzMREhKCHTt24Fe/+lWb+21ubkZzc7P9a7PZHNgnQuQH4pBDoClVuCr14rSkfEmbY/URwLWngXg3o+lyXvjUODVdjYFNLsH83MmZpobVPKmtrQUAxMfHOx2Pj4+331ZbW4s4cSrvFZ06dUJMTIz9nNYKCwsRFRVl/5OcnByA1hNpj5IbgHbk4nQiwvb/2PP+v29feaoT8/fUdKnEwOaubs0Ag65mOzoK5udOznQTjgLl8ccfR0NDg/3PkSNHlG4SkeKULlz1dhHzpL6r7f+ONUeAchc+d3ViSaYkRYrn1RjY5BLMz52c6SYcJSQkAADqWq18W1dXZ78tISEBJ1rtqXT58mWcPn3afk5rYWFhMJlMTn+Igp3ShatSLmLu2HuOHMKR0he+nNQcVOVXoSyvDOty1qEsrwyV+ZWKFc+rLbDJKZifO12lm5qjfv36ISEhAZ999hmGDx8OwFYftGPHDsyZMwcAMHr0aJw9exY7d+7EyJEjAQCff/45rFYrRo0apVjbibRGDYWrnhYDvH/E/S7rjYCr4SjeIRz5cwHB9pKrTkwquQv71SSYnzvZaCocNTY24scff7R/XVlZid27dyMmJgZ9+vRBQUEBnn32WVxzzTXo168fnnrqKfTq1QuTJ08GAKSmpmLcuHF44IEHsHr1aly6dAnz5s3D9OnTOVONyAdqKVx1dxEDgLe+fctlobMYjnpd7Ix3f/U2ept688LnhtoCm5yC+bmTxlbILi8vR0ZGRpvjeXl5WLNmDQRBwJIlS/Dmm2/i7NmzuPXWW/H666/j2muvtZ97+vRpzJs3D//4xz8QEhKCKVOm4NVXX0W3bt0ktYErZBN5X41bqZWdHYkF4wCc2njHT8CnfwHM/ZNg+j/WEBIFC1+u35oKR2rAcERk4y58qGmbBVdLDWQ2xqP0pTqgRw/g5EkFW0dEcmI4CiCGI6KrXIWPZFOy4vU7jtosUhl2LYy9egMGA9DSAnTSVHUBEbUTw1EAMRwROVNihewOuXwZCA0FBAGorQVarY1GRPrky/WbH5mIyImvYUdzhaudOgExMcCpU0B9PcMRaY7mPpBoEMMREdkptR2I7OLibOGorg4YOlTp1hBJFjSvUYXpZhFIIuoYJbcDkZ24jVB9vbLtIPJBUL1GFcZwRESKbwciF4vVgvKqclSHXQQAWOtc76lIpDbB8hpVC4YjIlJ8OxA5lBwoQUpxCjLWZuAfZ3YAAFb+czE/bZMmBMNrVE0YjohIFduBBFLr4Qhxlezw0+dcDkeIPUzr96xHeVU5P42T4vT+GlUbFmQTkWq2AwkEV8MRrTefLdhagEmDJsEYYmTBK6mSnl+jasRwRERI65OGJFOS1+1AxH3L2kOp6ceuhiPqu9r+H9vkPBxx+sJp5G7IbfMzEAte1bDqNwUnOV6jjoJ9uQCGIyKCMcSI4nHFyN2QCwMMLrcDKRpX1O43R1e9MT279sRvrv8NJg2aFNA3XlfDDGLPUVzT1WM15hos+myR24JXAwxOPUyB5O3CFOwXrmAU6NeoI/aesuaIiK7ISc3Bpmmb0NvU2+l4kimpQz0m7qYfnzx/EkX/LkLG2gykFKcErDDa1TCDfVjt/NVj9efrVVHw6lg4fnfJ3W1+Pt5uJ/0K1GvUEZcLsOH2IT7i9iGkd/7slbBYLUgpTvEYOoDAblYrtsFxOKL7eeD0ctvtYU8C8THJKLyjEL95/zde769gVAFeGfeKX9soEi9MrXuvxJ/PI2MewUvbX3J7O4f9gkOgeg69vV7FobvK/EpN9lRyb7UAYjgiLVDLsEt5VTky1mZIOjeQb7xi6ACuDJFZgZZngE4C0HsBsPL+vyEmPEZyW/827W8BC3GegqTRYIRFcD1zTusXLlKe1NdrWV6ZtrYMusKX6zeH1Yh0Rs5hF29T3n2ZVhzIYavWwxFCCFB/ZWjtL2NWICc1x17w6o1Ye+Rpen97lgLwto4NALfBCOA6N9RxXC7gKlXixhAAACAASURBVBZkE+mIu2GZQMy2klK02Z5pxYF6481JzcGkQZPsPWqR//MkcPAn3BF5PYCrBa9TNkzxeD+OIcTVp+f2FrP663kHw4WLAoPLBVzFniMinZBzewGpRZtib4xYEyNFIN94jSFGpKekY8awGejWu5/t4IkT9ttzUnNQcHOBpPtyFUI6Uszqr+cdDBcuJQTDwqDeXq8GGJBsSvbbcgFqxnBEpBNybS/gSwgTe2OkkP2NV9x81iEcAcCkQZMkfXvrENLRcColSBoNRl64FBAsMwQdX6+t/535e7kAtWM4ItIJueoFfA1hYr2Pp3oeRd54xXBUX+90uL2fnjsaTr1dmAwwYMHoBW5vB4LnwiWnYJvaLsdyAVrAcESkE3LVC7QnhOWk5qAqvwpleWUoGFWA2K6xTucq8sYrhqO6OqfD7f307I9w6u3CtPzO5bxwyUjOoWo1cXy9rstZh7K8MlTmVwbVvy8WZBNpmOOU/biIOCRFJqHmXGC3F2hvCBPrfdJT0vHS2JeUX2og9kpAa9VzBFwNKa4Kq4vGFbm8SPgrnLYuHG/98xFvL68qR3lVOQDYf67kX770Burt5y++XoMVwxGRRrmaFdUjvId9q4tAbS/gjz2eVPHG66bmSOQtpLTmz72vvP18Pjj0gdPv/tltzwbd9g5y4NT24MVhNSINclcHcfrCaQBATHiM03F/DrvopmjTTc2RI8fZbekp6R6fk1w/l2CrgVESp7YHL4YjIo3xVgdhgAHhncLx6b2fBqxeQBdFm25qjjoi0D+XYK2BUQqntgcvDqtpjFq2hSDlSKmDOHruKIwhRswYNiNg7fB12El1xJqj8+dtf7p29cvdBvLnEsw1MEoQewNzN+QGdKia1IfhSEPau/IuqU9HQq4SdRDu2quK2qH2iowEwsKA5mbb0Frfvn6760D9XFgDI7/2FOeT9jEcaYSc20JQYHU05Eqtb6hrqrMvxNgRug3lBoNtaO3IEdvQmh/DUaCwBkYZmu8lJZ8ZBEFoO3hNbvmyq6+/eNutm7txe6eW4Uh3IVfsopcScsV/D+5mRTnqaIjxR3tV7cYbgZ07gQ8/BLKzlW6NV95+91LfC9TyeiCSky/XbxZkq4SnfXvk2hZCr9Sy9L+/imk9zYpqrSMzmIKi+FesO3IznV9t/DEjTi2vByI1YzhSAW9vVqwzaD81TXv2Z8h1NyvK1X0C7QsxQRHKJUznV5uOzIhT0+uBSM0YjhQm5c2KdQbto7aeD3+HXHGJ/1eyXvF4XntDTFCE8gBM5/cHbzvAt2d7B7W9HojUjAXZCpKyXk3B1gL8+PCPflt5N5iobdpzIEKuMcSI+Ih4Sef6GmKCIpR72EJEKVIL4H2dEae21wORmrHnSEFS36y2H92ujxWJZeD4ifuznz6T9D1y9XwEakG5QIWYoFgAT2U1R4Ec9gqKnkAiP2E4UpAvb1a6WJE4wFrXbj277VlJ3ydXz0egtpcIVIjRzTYhnqio5ijQw15B0RNI5CcMRwry9c2qPXUGwcLdJ25PlOj5CETIDWSI0X0o97L5rJwCXQAfFD2BRH7CmiMFtWcXb02vSBwgnj5xu6Nkz0cgFpQL5Cq+ul4Az3FYTRBsC0MqJNDDXtwKg0g6hiMF8c3KP7x94nZF6aX/AxFyfQkxvi4CGKhQrvhihGI4amkBzp0DZFrY1RU5hr24FQaRNAxHCuObVcdJ/ST9ZNqTGBI7RF89H61ICTFq2Q5EFe2IiLD9aWqyTedXMBy1pye5PXTdE0jkJ9w+xEeB2j5E8U/QGlZeVY6MtRlezyvLKwv6IUm1bAeilnYAAPr3Byorgf/9X2DMGHke0w3x5wLAZU+yLuq8qF14jeg4X67fDEc+UmJvNfLMX/tN6Z1a9uhTSzvsRo0CvvoKeP99YPLkwD+eF6561JJNyexJDmKq6GXVAe6tRkElKKac+4FatgNRSzvsVDSdH+CsVHLGLV+UwXBEuqD7Ked+oJZFANXSDjsVbiEi1o7NGDYD6SnpQR/sgxW3fFEOC7JJN/RUaBqI+gK1LAKolnbYqazniEjELV+Uw3BEuiJ+4hbDxYZ9GzQXkgJVXyDXbCittMNOhfurEQEq7GUNIhxWI91pvY1IxtoMpBSnaGJsPpD1BWqpzVJLO+xUOKxGBKiwlzWIMByRrmi5eFGO+gK11GappR0AOKxGqsUtX5TDqfw+4lR+9VLdFHEfyblek1rWTFFFO3btAkaMAOLjgdpaeR+byAuufeU/vly/WXNEuqH14kU56wvUskefL+0IWJASe45OngSsViCEHeqkHtxFQRkMR6QbWi9e1Gp9gRy9PwFdBE8syLZYgDNngB49OnZ/FDTk6vnU00xcrWA4It3QargQqW4WlwRyrNzrbqsRsY6sw8MKoaFAdDRw9qyt7ojhiCSQe9VqtfT2Bgv2H5NuaL14UXWzuLyQo/hdtkXwxN6jEyc8tqW8qhzr96xHeVU5F94LYlqe+EHSMByRbmgtXLiiqllcHsgVWmTbakSsO3ITjrS8PAT5F1etDg4MR6QrWgkXnmhhby25QotsdWQepvOzl4AcqW5vQAoI1hyR7uiheFHt9QVyhRbZ6sjc9Bx56yUwwICCrQWYNGiSpv59UftpfeIHScNwRLqk9nChdXKFFtmK1N1sIaL15SHI/7Q+8YOk4bAaEflMruJ32erI3Gwhwl4Cak2NEz84WcD/GI6IyGdyFr/LUkfmpuaIvQTUWnv+7QcyvHCyQGBw+xAfcfsQoqtcrfWSbEoOyMq9AV1wr6wM+OUvgdRUYP9+p8dMKU5xO6wHADHhMdiQuwHpKemsOwoiUv/tB3I9JHdrgPl9axELgKYrfxod/t8I4HYA4R1/CDn4cv1mOPIRwxGRM1Xsj9ZRe/cCw4YBPXu26T1yt7dVa4FcAFArdPFvwQfenm8gw4urvSRDrCGIaIlAt5ZuiGyJRP/Q/vjwPz+E8YLROdQ4hpvWX7v6+0UPDTkM4Jp2PQXZMRwFEMMRkQ6dOGHbeNZgAFpagE7Oc1VcffpvLdg3ApV7xWi187QRdog1BN1auuHasGvx7+n/vhpe3IUTF6HGfNqM6ppqdGvphohLtkAUfjnAXTghALoBiHD4/3sArg3sw/oLw1EAMRwR6ZDFAnTuDAgCUFtrC0qtT7lSNzJt0zScvnDa5d2Is+cq8yt13WPSmmzDO3K7jLZhxVsvy5W/nzxxEt/99J0tvFzpzZErxFgMFjSGNqIxtBGRMZEwxZiuBhox1Dh+3Q2ub49sdTwMcFOHrgm+XL85lZ+IyGi07al28qRtWM1FODKGGGEMMboNRoD0qf16Gn5SxVpQlyFteMiHcINGAM3tb1JP9MQduMPjORaDBZYIC0IjQ9sGEU9hJgLY07gH+f+bj6bQJjSGNqKpc5M9EDV3araHmLK8Mi4z0Q4MR0REgG3G2smTtun8Q4e6PMUfU/v1Nvzk01pQSenSQonU4SU/hBhJjGg7nOSplyUCONR8CEu/WeoUXFoHmeZOzSib2b7wMsQ6BD80/aCpjaq1hOGIiAiwhaP9+11uISLq6NR+d8NP4lYkqhl+ugRpoaQR6FnZE698/4rTsJHjUJJ4PPq5aNv9BpIYYlwFFzc9MJLCTih8Hk4aaB2IfxX/y2N46ch6SOKSArkbcmGAwekxtLKXpJoFbTh67bXXsGLFCtTW1uKGG27AypUr8Ytf/ELpZhGRUrxsPgt0bMXugAw/tQ4xvva4uPueFmkPDwBDr/znk07w2OvitRbG3dftCDGBIkd4EdcAc9UTGYjlNIJJUIaj9957DwsWLMDq1asxatQoFBUVISsrC4cOHUKc+AZJRMHFzRYijtp9wWsBKvZXAEeAwS2D2xToOvayHP3xKPp26us91PgYYtqlE1wHF4f/WyOseP3A66gVatHYuRHnws45DSM1hTahW0w3lM4phTHSaPvesAC3WyXkCC962EtSjYJyttqoUaNw00034c9//jMAwGq1Ijk5GQ8//DAWLVrk8Xs5W41Ip55+GliyBLjvAeDZN2ExW7D7x904e+os4oQ4XBd+HUIuhACNwL7KfSjbWwZLo8UeanoKPXF9xPWItca2DTGBHk4Se2LcDSl5CDceh6BCpT28u7WgND9bzU/0VICvZZzK70FLSwu6du2KTZs2YfLkyfbjeXl5OHv2LD744AOn85ubm9HcfLXaz2w2Izk5meGISEktaN9idp56YsyrAescAJMBvB+QZjcbm9HU2dab0tS5qU0vS2NoI7JuyEJSryTPdTCtv5YYYgJJztXSidqDU/k9OHnyJCwWC+JbTdWNj4/HwYMH25xfWFiIZcuWydU8Iv0Q0DbEuBse8jXcXA5Eg68Mq+EEWkJaXAaXxtBGDB84HMmJyZJmLDn25FjCLRj4+kCv9Uoz82dqctdLpYZ32CtDgRB04chXjz/+OBYsWGD/Wuw5ItINxxDTnh4XT38PSIhxEArJU6u9Di/tiwNygcqeX6H/PNdFMR1Z5NEI/c8uMoYYZV1TR2/LIpB6BF046tmzJ4xGI+rq6pyO19XVISEhoc35YWFhCAsLkupBUjcBtvVcAtET479Nwl0Lg7ThIXezldyd49d3MFvPUfdz7hOd1EUe3eHsIv/RzLIIpElBF45CQ0MxcuRIfPbZZ/aaI6vVis8++wzz5s1TuHWkC44hxl9Tq8W/WwPc9jBIn1rtS7jRwjvNlZmq0c1A6GWgxUObpS4G6QpnF3WcKlblJl3TwluW3y1YsAB5eXm48cYb8Ytf/AJFRUVoamrCrFmzlG4ayUlKiJEabuQOMV3Qvv2RvK0ZE5TvCFdER0MwhsBgsSK2CaiJcn+q1MUg3ZF7+ElvfFqVmz9naoegfCv89a9/jfr6eixevBi1tbUYPnw4tm7d2qZIm1RCAHAR/h1GEvdYCvRwUhe4Dy6+TK1u3VsTlK/cAAsJAWLjgNpaxLsJR9ySQR38sY2LlrDoXH5B+xY7b948DqP5mwDgAjwHlPaGm0AvOCH2xHR0pd7W38/3L00xxNnCUex56LZoWg86uo2LlrDoXBlBG46Cmhhi/N0TI0eICYf3oOKtl6b18a5giCGbK3VHzwx9GPtC32fRtEp1ZBsXLWHRuXIYjvTuXgB70Lb3Rq4Q057aF3e3syeGAu3KFiI3dU5BVX4VhzJUKhg2XWXRubIYjvTuMIDvPNzeFT4vZue1XoY9MaRVDpvPsmha3fS+LAKLzpXFcKR3L8E2hOZqmKkrNLkSL1HAiOGo1TpopE56XhYh2IrO1YbhSO+0PeROV3C2ikzEcFRfr2w7SDK99vAFU9G5GvktHB04cADZ2dn46aef/HWXRATOVpHVlZojnDihbDso6AVL0bla+W1QpaWlBT///LO/7o6IcHW2SuvaA3G2SsmBEoVaplPsOSKVEIvOgatF5iK9FJ2rmeSeI8fNV12p55sJkV9xtooCWHNEKqL3onM1MwiCIGlSt9FoxPDhw2EymVze3tjYiG+//RYWS6CXHFaW2WxGVFQUGhoa3P4siPyhvKocGWszvJ5Xllemy5oLRZw7B4iv68ZGICJC2fYQgTWH/uLL9Vtyz9HAgQMxf/58/OY3v3F5++7duzFy5EjfWkpEbnG2igK6dQPCwoDmZlvdUb9+SreISLdF52omueboxhtvxM6dO93ebjAYILETiogk4GwVBRgMrDsiIuk9Ry+//DKam5vd3n7DDTfAag30VuRE/qfWLmvOVlFIXBxw5AhnrOmQWl/rpD6Sw1FCQgIAoKysDBkZrusg3njjDTz00EP+aRmRDNQ8TT4YtkhQJU7n1yU1v9ZJfXyeyj9u3Dg8+uijuHTpkv3YyZMnMXHiRCxatMivjSN9s1gtKK8qx/o961FeVQ6LVd5ifi1Mkxdnq/Q29XY6nmRK4qaTgcJhNd3Rwmtd75R+v/eVz4tAlpWV4b777kNpaSnWrVuHyspKzJ49G4MGDcLu3bsD0UbSIaU/xWlpmryet0hQJYf91Uj7tPRa1yul3+/bw+eeozFjxmD37t0YOnQoRowYgV/96leYP38+ysvL0bdv30C0kXRGDZ/ifNnUUQ3E2Sozhs1Aekq6xzdxpT+hKf34HcaeI13R2mtdb9Twft8e7do+5PDhw/jmm2+QlJSEY8eO4dChQzh//jwiuCYIeaGWT3F6nSav9Cc0pR/fL1hzpCu+vtZZtO0/anm/bw+fe45eeOEFjB49GnfeeSf27t2Lr776Crt27cL111+PioqKQLSRdEQtn+L0OE1e6U9oSj++37DnKODk7F305bVecqAEKcUpyFibgbtL7kbG2gykFKdo59+uyqjl/b49fA5HxcXF2Lx5M1auXIkuXbpg6NCh+Oqrr5CTk4P09PQANJH0RC09NuI0+dZ7FokMMCDZlKyZafLePqEBQMHWgoBdhLw9vgAB+VvztTHExi1EAkruACL1tV7fVK+PcK8ianm/bw+fw9GePXswfvx4p2OdO3fGihUr8Mknn/itYaRPaumx0dumjkp/QvP2+ABw1HwUz217LiCP71fisFp9PcCFbf3K195Ff/QwSXmt/9fY/8KCTxYo9uFCr9Tyft8ePoejnj17ur3t9ttv71BjSP/U1GOjp2nySn9Ck3q/S8qXqP8TuBiOWlpse62RX/jau+nPHiZvr/WeET01O/yjZmp6v/dVuwqyidpLbQsb6mWavNKf0Hy5X7UWYNp17WrbY62x0Ta0xg2m/cKX3s3TF04jd0NumyAl9jC158OLp9f6+j3rJd2HGod/1Ext7/e+8LnniKij1NZj48s0ebVS+hOa+PhSaOITOIuy/U5qsKgx1wSsfs7da13pDxd6prb3e6nYc0SK0EuPjVo4fkJrTY5PaOLjT9kwRdL5qv8EHhsL/PQTp/P7kdRgUX++XnIPk792quc+hoGlxfd7hiNSjPgpjvwnJjwGpy6canPszYlvBvwTWk5qDpalL8OS8iVez1X9J3D2HPmd1AAS2zVW0v35M2BrefhHK7T2fs9hNSIJ1L7qszgLqHUwAuDyWKD8Ke1P6B3Z2+3tai7AdMLp/JL48rqQOkO09fCLO/4O2Fod/qHAYM8RkRdqX/XZ0ywgALKuQmsMMeLV8a/ah/c0+wmcPUdeted1IQYQV99XNK4IOak5sFgtig1xaXH4hwLDIAhcyMMXZrMZUVFRaGhogImzWHRP7JFp/SYtXujV8ImyvKocGWszvJ5XllcmW7e2qwtnsinZfgFUvVdeARYsAKZPB9ZLm8mkRoHaCqOjrwtv7RLvH3AdsNXwuiPt8eX6zZ4jIje0si+Q0mscuaL5T+A66DkKVI+nP14X3upPpPQwEQUSwxGRG76sy6JkoaFapyFrrQDTiRiONDpbzV3PzlHz0XavEySS63Wh+YBNmsZwROSGGntkXOE05ADQcDjyVoMmQOhQj2d7XhftHd7TdMAmTeNsNSI31Noj05qUWUD3j7gfG/ZtUOVMO1UStxA5eRKwWpVti4+k7HPXkYU4fX1dcKd70iKGIyI3lF512hfupiHHhMcgJjwGS8qX8MLkCzEcWSzA6dPKtsVHNeYav57Xmi+vC183miVSC4YjIjekrsuilhqInNQcVOVXoSyvDOty1mFZ+jKcunCqzTpHvDBJ0Lkz0L277e8aK8quPy+tvVLPa03q6wJAwLYBIQo0hiMiD7S2MJxYozHtuml469u3XJ7DC5NEYu+RxuqOpK4wLfU8V6S8Lnwp3CZSGxZkE3mhxVkzWplpp2pxccDhw5rrOZK6wrTU89zx9rrQyoQGIlcYjogk0NqsGV6Y/EDsOdLYFiJiTZCncOyvWjlPrwutTGgIVoFaIFQvGI6IdIgXJj/Q6EKQjpuoAspt4cIlJtRL7VsiqQFrjoh0SEsz7VRL7DnSWDgC1FErp7UJDcGCMwilYTgi0iFemPxA7DnS2LCaqPXsxbK8MlTmV7YJRharBeVV5Vi/Z73f18FSQ0ijq7xt/QJwooaIw2pEOsX9qTpIo8NqjrzVyskxvKLFCQ16xYka0jEcEbWip0JFXpg6QMNbiEjhbv81cXjFnz07WpvQoFecqCEdwxGRAz0WKnq7MOkpDPqVhmuOvJE6vNLe/ddInThRQzqGI6Ir5PwkrRab9m3C77f83mm1ZDWHQVmDnNhzdOoUcOmSbdVsnfBl/zX2+OgHZxBKx4JsIgRnoeJjpY9h6qapbbaROGo+qspZK7JvYNqjB2C4Usx+6pTnczVG6rDJBwc/CHBLSE6cqCEdwxERfCtU1ION+zZixfYVbm8XIKgqDCoy/dhotAUkQHd1R1KHTf6656+q+TdA/sEZhNJwWI0IwVWoaLFa8Pstv/d6nhLDKq6GzQDPG5gaYAhcfUxcHHDypO7qjtL6pKFn1544ef6kx/Pqz9dzaE2HOFHDO4YjIgRXoeK26m1eL4oiOcOgu2L4B0Y8oNz047g4YP9+3fUcGUOM+M31v0HRv4u8nquHDwS+CJYJCpxB6BnDERGCq1DRl4udXGHQUzH8kvIlku4jIBdxjS8E6cmkQZMkhSM9fCCQSo+zVal9WHNEhOAqVJR6sYvtGitLGJRSDC9FQC7iOp7Ozy1mnHFbDXLEcER0RbAUKooXRW9ev+t1WcKglGnlngT0Iq7jhSCD6QOBN8E4W5U8YzgiciB1PyotEy+K7noMAODRMY8i97pcWdrjy3CY7BdxHWwh4kmwfCDwJthmq5J3rDkiaiUYChXd7bsW2zUWr931GqZeN1W2tkgdDluWvgxvffuWvPvE6TwcAZy5BEhfzynYitODmUEQBOmD+gSz2YyoqCg0NDTAZDIp3RyiDlHDzByL1YKU4hSvxfCV+ZUAIG97t20DbrsNuOYa4PDhwD0OKcZitSDx5cQ2i6G6UpZXpvsPTnrmy/WbPUdEQUwNvWTiMF/uhlwYYHAKSK6GzWRtr45nq5HNtuptkoKRXBMUSB1Yc0REilNt7YsYjsxmoLlZmTZQQEkdKrvn+nuCaqgx2LHniIhUQZW1L9HRQKdOwOXLtrqjJO+z/EhbpNa8TRo0KcAtITVhOCIi1VDDMJ8Tg8G21tHx47bp/AxHuhOIBWDVUMtHHcNwRETkSVycLRzpeMZaMPO15s0brrKtD6w5IiLyhEXZuuevmjeusq0f7DkiIvJEx1uI0FUdrXnztsq2AQYUbC3ApEGTOMSmAQxHRK2wXoCc6HgLEXLWkZo3X1bZVlVdHbnEcETkgPUC1EYQrJJNHSd1SQCusq0Nmqk5eu655zBmzBh07doV0dHRLs+prq5GdnY2unbtiri4ODz66KO4fPmy0znl5eUYMWIEwsLCMHDgQKxZs0aG1pMWsF6AXGLPEUkgdUkAqeeRsjQTjlpaWjB16lTMmTPH5e0WiwXZ2dloaWnB9u3bsXbtWqxZswaLFy+2n1NZWYns7GxkZGRg9+7dKCgowP3334+PP/5YrqdBKsVduckt9hyRBOKSAO42dDbAgGRTMlfZ1gjNhKNly5Zh/vz5GDZsmMvbP/nkE+zfvx/vvvsuhg8fjvHjx+OZZ57Ba6+9hpaWFgDA6tWr0a9fP7z88stITU3FvHnzkJubi1deeUXOp0IqxF25yS2xIJs9R+SBuCQAgDYBqT1LApCyNBOOvKmoqMCwYcMQHx9vP5aVlQWz2Yx9+/bZz8nMzHT6vqysLFRUVMjaVlIf1guQW5zKTxKpdhsc8pluCrJra2udghEA+9e1tbUezzGbzbhw4QLCw8Pb3G9zczOaHfZUMpvN/m46qQDrBcgtMRxduAA0NQEREcq2h1RNldvgkM8U7TlatGgRDAaDxz8HDx5UsokoLCxEVFSU/U9ycrKi7aHAYL0AuRURAXTpYvs7h9ZIAnFJgBnDZiA9JZ3BSIMU7TlauHAhZs6c6fGc/v37S7qvhIQEfPXVV07H6q50gyckJNj/X9eqa7yurg4mk8llrxEAPP7441iwYIH9a7PZzICkQ/7eQoC0x+36VgaDrfeoutpWlN2vn9JNJaIAUzQcxcbGIlYsduyg0aNH47nnnsOJEycQd6UbvLS0FCaTCUOGDLGfs2XLFqfvKy0txejRo93eb1hYGMLCwvzSRlI3sV7A1TpHReOKWC+gY17XtxLDkcbqjrigKVH7aKbmqLq6GqdPn0Z1dTUsFgt2794NABg4cCC6deuGsWPHYsiQIbj33nuxfPly1NbW4sknn8TcuXPt4eZ3v/sd/vznP+Oxxx7Db3/7W3z++efYsGED/vnPfyr51EhFWC8QfMT1rVov4yCub7Vp2ibkaHALES5oStR+BkEQ2i7sokIzZ87E2rVr2xwvKytDeno6AODnn3/GnDlzUF5ejoiICOTl5eGFF15Ap05XM2B5eTnmz5+P/fv3IykpCU899ZTXoT1HZrMZUVFRaGhogMlk6ujTIiIFWawWpBSnuF3GwQADkkxJqNqVjpB33gEKC4FFi2Rupe/cBT5xiJgzpygY+XL91kw4UguGIyL9KK8qR8baDK/n/Vzza/R56z1gwQLg5ZdlaFn7SQ18lfmV7BEl1QnkULAv12/NDKsREfmb1HWrTnULQR9AE7PVuAEqaZWahoJ1swgkEZGvpK5bFd6rj+0vGghHXNCUtEhte1syHBFR0JK6vtU1g2+xHdBAQTYXNCWtUePelgxHRCpnsVpQXlWO9XvWo7yqnJvf+pHjfliuCBAwfeh0GBOuBAkNTOXngqakNWrc25LhiEjFSg6UIKU4BRlrM3B3yd3IWJuBlOIU2buY9SwnNQePjHnE7e0vbX8JHzXstH1RXw+ofA4LN0AlrVHjUDDDEZFKqW0MXq8sVgvW713v8Zw/fPOM7S+XLgENDTK0qmO4ASppiRqHgjmV30ecyk9y4HRs+Uidzn9pRTg6NV0ADh8GrrlGhpZ1HFfIJi0Q3+9qzDUu64789X7ny/WbPUdEM4WLpwAAFdtJREFUKqTGMXi9ktpVf6F7pO0vGqg7EnEDVNICNQ4FMxwRqZAax+D1SmpXvTW2p+0vGpixRqQ1ahsK5iKQRCqkxjF4vRJnd3nr0jf17g/s2s9wRBQgatrbkj1HRCrE6djykdqlb4iPtx3U0LAakdaoZSiY4YhIhdQ4Bq9nkrr04+JsB9lzRKR7HFYjUinxgu1qr6GicUWcju1nXrv0Y2Nt/9fAFiJE1DEMR0QqpqYx+GAgdum7xJ4joqDBcESkch4v2CQfMRyx5ohI91hzREQkBXuOiIIGwxERkRRizdHJk4DVqmxbiCigGI6IiKQQw5HVCpw+rWxbiCigGI6IiKTo3Bno3t32dw6tEekawxERkVRi3RGn8xPpGsMREZFU4tAae46IdI1T+YlI1SxWi3rWeeJ0fqKgwHBEFCCquqhrVMmBEpcrhBePK1ZmhXBO5ycKCgxHRAGguou6BpUcKEHuhlwIEJyO15hrkLsh9+qeZ3LiFiJEQYE1R0R+Jl7UHYMRcPWiXnKgRKGWaYfFakH+1vw2wQiA/VjB1gJYrBZ5G8aeI6KgwHBE5EeqvahrzLbqbW3CpSMBAo6Yj2Bb9TYZWwXOViMKEgxHRH6k2ou6xhw/d9yv5/kNwxFRUGA4IvIj1V7UNSYxMtGv5/kNp/ITBQUWZBP5kWov6hqT1icNSaYk1JhrXA5RGmBAkikJaX3S5G2Y2HN06hRw6ZJt1WyN4SxKIu8YjkhXlH7jV+1FXWOMIUYUjytG7oZcGGBw+lkaYAAAFI0rkv+iHhMDhITY9lc7dQpISJD38TuIsyiJpOGwGulGyYESpBSnIGNtBu4uuRsZazOQUpwi6+ww8aIOXL2IixS9qGtQTmoONk3bhN6m3k7Hk0xJykzjBwCjEejRw/Z3jdUdcRYlkXQGQRDafrwlt8xmM6KiotDQ0ACTyaR0c+gKd2viiIFE7oupq0/oyaZkFI0r4id0HyndG9jG/fcDFy4AS5YA116rXDt8YLFakFKc4naygNijWZlfyeBOuuXL9ZvhyEcMR+qj1jd+1V3UKWiVV5UjY22G1/PK8sqQnpIe+AYRKcCX6zdrjkjzfJk+L+cbvzHEyAsNqYLSsyj5QYG0huGINE/pN34itVNyFiWLwEmLWJBNmsfp80SeibMoW08SEBlgQLIp2e+zKFkETlrFcESap9QbP5FWKDGLklvpkJYxHJHmcfo8kXdyL43ArXRIy1hzRLogvvG7qm3g9Hkim5zUHEwaNEmW4mjWApKWMRyRbsj5xk+kVXLNomQtIGkZwxHpCqfPE6kDt9IhLWPNERER+R1rAUnLGI6IiCggVLk/HpEE3D7ER9w+hIjIN1whm9SA24cQEZFqsBaQtIbDakREREQOGI6IiIiIHDAcERERETlgOCIiIiJywHBERERE5IDhiIiIiMgBwxERERGRA4YjIiIiIgcMR0REREQOGI6IiIiIHDAcERERETlgOCIiIiJywHBERERE5IDhiIiIiMgBwxERERGRA4YjIiIiIgcMR0REREQOGI6IiIiIHDAcERERETnQRDiqqqrC7Nmz0a9fP4SHh2PAgAFYsmQJWlpanM77/vvvkZaWhi5duiA5ORnLly9vc18bN27E4MGD0aVLFwwbNgxbtmyR62kQERGRBmgiHB08eBBWqxVvvPEG9u3bh1deeQWrV6/GE088YT/HbDZj7Nix6Nu3L3bu3IkVK1Zg6dKlePPNN+3nbN++HTNmzMDs2bOxa9cuTJ48GZMnT8bevXuVeFpERESkQgZBEASlG9EeK1aswKpVq/DTTz8BAFatWoU//elPqK2tRWhoKABg0aJF2Lx5Mw4ePAgA+PWvf42mpiZ8+OGH9vu5+eabMXz4cKxevVrS45rNZkRFRaGhoQEmk8nPz4qIiIgCwZfrtyZ6jlxpaGhATEyM/euKigrcdttt9mAEAFlZWTh06BDOnDljPyczM9PpfrKyslBRUeH2cZqbm2E2m53+EBERkX5pMhz9+OOPWLlyJR566CH7sdraWsTHxzudJ35dW1vr8RzxdlcKCwsRFRVl/5OcnOyvp0FEREQqpGg4WrRoEQwGg8c/4pCYqKamBuPGjcPUqVPxwAMPBLyNjz/+OBoaGux/jhw5EvDHJCIiIuV0UvLBFy5ciJkzZ3o8p3///va/Hzt2DBkZGRgzZoxToTUAJCQkoK6uzumY+HVCQoLHc8TbXQkLC0NYWJjX50JERET6oGg4io2NRWxsrKRza2pqkJGRgZEjR+Ltt99GSIhzp9fo0aPxpz/9CZcuXULnzp0BAKWlpRg0aBC6d+9uP+ezzz5DQUGB/ftKS0sxevRoPz0jIiIi0jpN1BzV1NQgPT0dffr0wUsvvYT6+nrU1tY61QrdfffdCA0NxezZs7Fv3z689957KC4uxoIFC+zn5OfnY+vWrXj55Zdx8OBBLF26FN988w3mzZunxNMiIiIiFVK050iq0tJS/Pjjj/jxxx+RlJTkdJu4EkFUVBQ++eQTzJ07FyNHjkTPnj2xePFiPPjgg/Zzx4wZg3Xr1uHJJ5/EE088gWuuuQabN2/G0KFDZX0+REREpF6aXedIKVzniIiISHuCYp0jIiIiokBgOCIiIiJywHBERERE5IDhiIiIiMgBwxERERGRA4YjIiIiIgcMR0REREQOGI6IiIiIHDAcERERETlgOCIiIiJyoIm91YiI9MxitWBb9TYcP3cciZGJSOuTBmOIUelmEQUthiMiIgWVHChB/tZ8HDUftR9LMiWheFwxclJzFGwZUfDisBoRkUJKDpQgd0OuUzACgBpzDXI35KLkQIlCLSMKbgxHREQKsFgtyN+aDwFCm9vEYwVbC2CxWuRuGlHQYzgiIlLAtuptbXqMHAkQcMR8BNuqt8nYKiICGI6IiBRx/Nxxv55HRP7DcEREpIDEyES/nkdE/sNwRESkgLQ+aUgyJcEAg8vbDTAg2ZSMtD5pMreMiBiOiIgUYAwxonhcMQC0CUji10XjirjeEZECGI6IiBSSk5qDTdM2obept9PxJFMSNk3bxHWOiBRiEASh7TxScstsNiMqKgoNDQ0wmUxKN4eIdIArZBMFni/Xb66QTUSkMGOIEekp6Uo3g4iu4LAaERERkQOGIyIiIiIHHFYjVWHtBRERKY3hiFSDu5MTEZEacFiNVIG7kxMRkVowHJHiuDs5ERGpCcMRKY67kxMRkZowHJHiuDs5ERGpCcMRKY67kxMRkZowHJHiuDs5ERGpCcMRKY67kxMRkZowHJEqcHdyIiJSC4MgCG3nT5NbvuzqS77jCtlERBQIvly/uUI2qQp3JyciIqVxWI2IiIjIAcMRERERkQOGIyIiIiIHDEdEREREDhiOiIiIiBwwHBERERE5YDgiIiIicsBwREREROSA4YiIiIjIAVfI9pG424rZbFa4JURERCSVeN2Wsmsaw5GPzp07BwBITk5WuCVERETkq3PnziEqKsrjOdx41kdWqxXHjh1DZGQkDAaD0s3xymw2Izk5GUeOHOFGuSrC34s68feiTvy9qJPWfi+CIODcuXPo1asXQkI8VxWx58hHISEhSEpKUroZPjOZTJr4xxts+HtRJ/5e1Im/F3XS0u/FW4+RiAXZRERERA4YjoiIiIgcGJcuXbpU6UZQYBmNRqSnp6NTJ46iqgl/L+rE34s68feiTnr9vbAgm4iIiMgBh9WIiIiIHDAcERERETlgOCIiIiJywHBERERE5IDhSAeqqqowe/Zs9OvXD+Hh4RgwYACWLFmClpYWp/O+//57pKWloUuXLkhOTsby5cvb3NfGjRsxePBgdOnSBcOGDcOWLVvkehq69Nxzz2HMmDHo2rUroqOjXZ5TXV2N7OxsdO3aFXFxcXj00Udx+fJlp3PKy8sxYsQIhIWFYeDAgVizZo0MrQ8ur732GlJSUtClSxeMGjUKX331ldJN0rUvv/wSEydORK9evWAwGLB582an2wVBwOLFi5GYmIjw8HBkZmbihx9+cDrn9OnTuOeee2AymRAdHY3Zs2ejsbFRzqehK4WFhbjpppsQGRmJuLg4TJ48GYcOHXI65+LFi5g7dy569OiBbt26YcqUKairq3M6R8p7mtoxHOnAwYMHYbVa8cYbb2Dfvn145ZVXsHr1ajzxxBP2c8xmM8aOHYu+ffti586dWLFiBZYuXYo333zTfs727dsxY8YMzJ49G7t27cLkyZMxefJk7N27V4mnpQstLS2YOnUq5syZ4/J2i8WC7OxstLS0YPv27Vi7di3WrFmDxYsX28+prKxEdnY2MjIysHv3bhQUFOD+++/Hxx9/LNfT0L333nsPCxYswJIlS/Dtt9/ihhtuQFZWFk6cOKF003SrqakJN9xwA1577TWXty9fvhyvvvoqVq9ejR07diAiIgJZWVm4ePGi/Zx77rkH+/btQ2lpKT788EN8+eWXePDBB+V6CrrzxRdfYO7cufj3v/+N0tJSXLp0CWPHjkVTU5P9nPnz5+Mf//gHNm7ciC+++ALHjh1DTk6O/XYp72maIJAuLV++XOjXr5/969dff13o3r270NzcbD/2xz/+URg0aJD962nTpgnZ2dlO9zNq1CjhoYceCnyDde7tt98WoqKi2hzfsmWLEBISItTW1tqPrVq1SjCZTPbf1WOPPSZcd911Tt/361//WsjKygpso4PIL37xC2Hu3Ln2ry0Wi9CrVy+hsLBQwVYFDwDC+++/b//aarUKCQkJwooVK+zHzp49K4SFhQnr168XBEEQ9u/fLwAQvv76a/s5H330kWAwGISamhr5Gq9jJ06cEAAIX3zxhSAItt9B586dhY0bN9rPOXDggABAqKioEARB2nuaFrDnSKcaGhoQExNj/7qiogK33XYbQkND7ceysrJw6NAhnDlzxn5OZmam0/1kZWWhoqJCnkYHoYqKCgwbNgzx8fH2Y1lZWTCbzdi3b5/9HP5eAqelpQU7d+50+hmHhIQgMzOTP2OFVFZWora21ul3EhUVhVGjRtl/JxUVFYiOjsaNN95oPyczMxMhISHYsWOH7G3Wo4aGBgD4/+3dXUgUbRQH8H+vuosmq4a6G4mgWEoZWUYxJaRYaXURXURJiRZamSKGCBVIdZUX0o100U16E0gFEWRZi7tLVCYoqymaUKmhtksfLAotqO15L8Jhpg+zt95dXf8/GHBnnhnOzJl55uz4DKveS7q6ujA9Pa3LS1paGhITE3V5+VWfthiwOApCr169QkNDA06ePKnOc7lcupMVgPrZ5XLN2WZ2Of19f5KXiYkJeL1e/wQaxD58+IAvX77w3F9AZo/7XDlxuVyIj4/XLQ8NDcWKFSuYt7/A5/OhqqoK27dvR3p6OoCvx9xgMHw3fvLbvPyqT1sMWBwtYGfPnsWyZcvmnF6+fKlbZ2xsDPn5+Th48CBKS0sDFHlw+y95ISJaTMrLy9HX14fm5uZAhxIQwfVjKEGmuroaxcXFc7ZJTk5W/x4fH0dOTg62bdumG2gNABaL5bs3CmY/WyyWOdvMLqevfjcvc7FYLN+9FTXfvJhMJoSHh88zavqZ2NhYhISE8NxfQGaPu9vtxsqVK9X5brcbGRkZaptvB8zPzMzg06dPzNsfqqioUAe4JyQkqPMtFgumpqbg8Xh0T4+018p8+rTFgE+OFrC4uDikpaXNOc2OIRobG0N2djYyMzPR2NiIf/7Rp1ZRFDx+/BjT09PqPKvVitTUVMTExKht2tradOtZrVYoivI/7+ni8jt5+RVFUdDb26vr5K1WK0wmE9auXau2YV7+PwaDAZmZmbpj7PP50NbWxmMcIElJSbBYLLqcTExMoKOjQ82JoijweDzo6upS29hsNvh8PmzdutXvMQcDEUFFRQXu3LkDm82GpKQk3fLMzEyEhYXp8jI4OIi3b9/q8vKrPm1RCPSIcPpzo6OjkpKSIrm5uTI6Oirv3r1Tp1kej0fMZrMUFhZKX1+fNDc3S0REhFy7dk1t8/TpUwkNDZX6+noZGBiQCxcuSFhYmPT29gZit4LCyMiIOJ1OuXTpkkRGRorT6RSn0ymTk5MiIjIzMyPp6emye/du6e7ultbWVomLi5Nz586p23jz5o1ERERITU2NDAwMyNWrVyUkJERaW1sDtVtBp7m5WYxGozQ1NUl/f7+cOHFCoqOjdW/c0N81OTmpXg8A5MqVK+J0OmVkZEREROrq6iQ6Olru3r0rL168kP3790tSUpJ4vV51G/n5+bJx40bp6OiQJ0+eyOrVq6WgoCBQu7TolZWVSVRUlDgcDt195PPnz2qbU6dOSWJiothsNuns7BRFUURRFHX5fPq0xYDFURBobGwUAD+ctHp6eiQrK0uMRqOsWrVK6urqvtvWzZs3Zc2aNWIwGGTdunXS0tLir90ISkVFRT/Mi91uV9sMDw/Lnj17JDw8XGJjY6W6ulqmp6d127Hb7ZKRkSEGg0GSk5OlsbHRvzuyBDQ0NEhiYqIYDAbZsmWLPH/+PNAhBTW73f7Da6OoqEhEvr7OX1tbK2azWYxGo+Tm5srg4KBuGx8/fpSCggKJjIwUk8kkx44dU7940O/72X1E2994vV45ffq0xMTESEREhBw4cED3RVxkfn3aQrdMRMSPD6qIiIiIFjSOOSIiIiLSYHFEREREpMHiiIiIiEiDxRERERGRBosjIiIiIg0WR0REREQaLI6IiIiINFgcEREREWmwOCIi+obD4cCmTZtgNBqRkpKCpqamQIdERH7E4oiISGNoaAj79u1DTk4Ouru7UVVVhZKSEjx8+DDQoRGRn/DnQ4hoSXn//j3Wr1+PyspKnD9/HgDw7NkzZGdn48GDB3j06BFaWlrQ19enrnP48GF4PB60trYGKmwi8iM+OSKiJSUuLg7Xr1/HxYsX0dnZicnJSRQWFqKiogK5ublob2/Hzp07devk5eWhvb09QBETkb+FBjoAIiJ/27t3L0pLS3HkyBFs3rwZy5cvx+XLlwEALpcLZrNZ195sNmNiYgJerxfh4eGBCJmI/IhPjohoSaqvr8fMzAxu3bqFGzduwGg0BjokIlogWBwR0ZL0+vVrjI+Pw+fzYXh4WJ1vsVjgdrt1bd1uN0wmE58aES0R/LcaES05U1NTOHr0KA4dOoTU1FSUlJSgt7cX8fHxUBQF9+/f17W3Wq1QFCVA0RKRv/FtNSJacmpqanD79m309PQgMjISO3bsQFRUFO7du4ehoSGkp6ejvLwcx48fh81mQ2VlJVpaWpCXlxfo0InID1gcEdGS4nA4sGvXLtjtdmRlZQEAhoeHsWHDBtTV1aGsrAwOhwNnzpxBf38/EhISUFtbi+Li4sAGTkR+w+KIiIiISIMDsomIiIg0WBwRERERabA4IiIiItJgcURERESkweKIiIiISIPFEREREZEGiyMiIiIiDRZHRERERBosjoiIiIg0WBwRERERabA4IiIiItJgcURERESk8S8WLupJjsKd8AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loss in Information due to\n",
      "dimensionality reduction: 62.484107971191406\n",
      "(high - we're losing real info)\n"
     ]
    }
   ],
   "source": [
    "# Let us now try a similar experiment with\n",
    "# uncorrelated data \n",
    "N = 100\n",
    "torch.manual_seed(42)\n",
    "x_0 = torch.normal(0, 100, (N,))\n",
    "x_1 = torch.normal(0, 100, (N,))\n",
    "X = torch.column_stack((x_0, x_1))\n",
    "\n",
    "# Perform PCA\n",
    "# Note that principal values are close to\n",
    "# each other - the spread is comparable in both\n",
    "# directions.\n",
    "principal_values, principal_components = pca(X)\n",
    "principal_values = principal_values.real\n",
    "principal_components = principal_components.real\n",
    "print(\"Principal values are {}\".format(principal_values))\n",
    "\n",
    "# Find the index with highest principal value\n",
    "major_index = torch.argmax(principal_values)\n",
    "minor_index = torch.argmin(principal_values)\n",
    "\n",
    "# Plot\n",
    "plt.figure()\n",
    "plt.scatter(X[:, 0], X[: , 1],\n",
    "            color=\"green\")\n",
    "plt.title('Uncorrelated dataset')\n",
    "plt.xlabel('x0')\n",
    "plt.ylabel('x1')\n",
    "\n",
    "# Let us plot the principal components\n",
    "draw_line(principal_components[:, major_index].numpy(),\n",
    "          min_x=-200, max_x=200)\n",
    "draw_line(principal_components[:, minor_index].numpy(),\n",
    "          min_x=-150, max_x=150, color=\"red\")\n",
    "plt.show()\n",
    "\n",
    "first_principal_vec = principal_components[:,\n",
    "                                           major_index]\n",
    "first_principal_vec = first_principal_vec.reshape((-1,\n",
    "                                                   1))\n",
    "X_proj = torch.matmul(X, first_principal_vec)\n",
    "\n",
    "# Information lost due to dimensionality reduction.\n",
    "# (here we've lost real info)\n",
    "X_back_proj = torch.matmul(X_proj,\n",
    "                     torch.linalg.pinv(first_principal_vec))\n",
    "info_loss = torch.sqrt(torch.mean((X_back_proj - X)**2))\n",
    "print(\"Loss in Information due to\\ndimensionality\"\n",
    "      \" reduction: {}\\n(high - we're \"\n",
    "      \"losing real info)\".format(info_loss))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
